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Abstract 

With applications in astroparticle physics in mind, we generalize a method for the solution 
of the nonlinear, space homogeneous Boltzmann equation with isotropic distribution function 
to arbitrary matrix elements. The method is based on the expansion of the matrix element 
in terms of two cosines of the "scattering angles". The scattering functions used by previous 
authors in particle physics for matrix elements in Fermi-approximation are retrieved as lowest 
order results in this expansion. The method is designed for the unified treatment of reactive 
mixtures of particles obeying different scattering laws, including the quantum statistical terms 
for blocking or stimulated emission, in possibly large networks of Boltzmann equations. Al- 
though our notation is the relativistic one, as it is used in astroparticle physics, the results can 
also be applied in the classical case. 
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1 Introduction 



Non-equilibrium processes in astroparticle physics, such as Big Bang nucleosynthesis (BBN), 
neutrino decoupling and more speculative ones, like baryogenesis through leptogenesis or the 
freeze-out of hypothetical relic particles, [1-7] are usually computed through the solution of the 
corresponding coupled system of Boltzmann equations [8-12] describing the time evolution of the 
one-particle distribution functions f l (t, k). Typically, in cosmology, it is anticipated that the rele- 
vant particle distributions are in exact kinetic equilibrium and of Maxwell-Boltzmann type. These 
assumptions, together with others, allow the Boltzmann equation to be linearized and integrated, 
which leads to coupled sets of chemical rate equations (mostly themselves dubbed Boltzmann 
equations in this context). This procedure drastically simplifies the numerical computation of the 
particle abundances, such that even the approximate solution of very large networks of Boltzmann 
equations, as in the case of BBN, becomes possible. However, in doing so, one looses the spec- 
tral information contained in the definition of the distribution functions and other fundamental 
properties of the Boltzmann equation are neglected as well. It is well-known that the solution of 
the full Boltzmann equations can lead to relevant corrections to the equilibrium results in several 
cases [15-17]. In the era of precision cosmology the inclusion of such non-equilibrium effects 
gains in importance. Regarding the use of classical kinetic theory for the description of phenom- 
ena in the (very) early universe there are concerns, originating in the belief that these calculations 
should be performed in the framework of non-equilibrium quantum field theory. These concerns 
are supported by recent results, revealing differences between the two approaches for simple toy 
models, at least in extreme non-equilibrium situations, see e.g. [13, 14]. However, it seems natural 
to attempt to include quantum effects in modified effective kinetic equations. Boltzmann equa- 
tions will continue to play an important role in cosmology at least at the relatively low energies of 
neutrino decoupling or nucleosynthesis, where the standard calculations give already quite good 
results. 

In general, a network of Boltzmann equations can be written as: 

L[f} = Y,C il [f 1 ,...,f,...,f N ], (1) 
i 

where there is one equation for each of the N participating particle species (i = 1 . . . N) and one 
collision term C %1 for each interaction with particles of the same and of other species. L denotes 
the Liouville operator, most commonly in Minkowski space-time, 

H/1<*) = ^, (2) 
or in Robertson-Walker space-time, 

mw= *JML- ai *£m, (3) 

with Hubble rate H = a /a. By writing the collision integrals as C ll [f l . . . f l . . . f N ], we have 
formally taken the possibility of multi-particle scattering processes into account. Usually only 
decays, inverse decays and 2 — 2 scattering processes, a + b <-> A + B, are considered. For the 
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latter ones the collision integral reads 



C al [f a f b f A f B ](k) 
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where we have used the short-hand notation f l k = f l {t, k) and to specify the quantum statistics 
of particle species i, i.e. £ l = +1 for Fermi-Dirac, £ l = — 1 for Bose-Einstein and = for 
Maxwell-Boltzmann statistics. \M\ 2 denotes the invariant matrix element squared and averaged 
over initial and final spin states. Note that we take 1 | 2 to include possible symmetrization factors 
of 1/2 for identical particles in the initial or final state. There is a vast number of different methods 
for the solution of the Boltzmann equation (mostly applied in different fields of physics), out of 
which only a few exploit the homogeneity and isotropy as imposed by the cosmological principle. 
So-called direct integration methods, where the collision term (4) is integrated numerically, seem 
to be most advantageous because they are characterised by high precision This is desirable since 
one wants to keep track of only small deviations from equilibrium. The direct numerical solution 
using deterministic methods is numerically expensive, mainly because of the multiple integrals 
in the collision terms. Therefore, it relies on the successful reduction of the collision integral for 
isotropic distribution functions. 

In the present paper a technique for this reduction of the collision integral is presented which 
generalizes previous results in high energy and astro physics [15, 18] for matrix elements in Fermi- 
approximation to (in principle) arbitrary matrix elements, relying on a series expansion of the 
matrix element. The resulting reduced Boltzmann equation contains only a two-fold integral over 
the magnitudes of the post-collisional momenta. The method is applicable to Boltzmann equations 
with and without quantum statistical terms and can be used independent of the dispersion relation, 
i.e. it can be used for massive and massless relativistic particles as well as for non-relativistic ones. 
The loss and gain terms can be treated collectively or independently. Thus the method represents 
an approach to treat reactive mixtures of all kinds of particles with different interactions, in a 
unified manner. 

The outline is as follows: In section 2 we show how the nine-dimensional collision integral for 
2 — 2 scattering processes can be reduced to a two-dimensional one, integrating out the energy 
and momentum conserving ^-functions. (Collision integrals for decays and inverse decays can 
be integrated in the same way.) In doing so, a certain angular integral over the matrix element 
arises. In section 3 we establish a simple numerical model for the reduced Boltzmann equation. 
The integral of section 2 is solved by expanding the matrix element in terms of the cosines of 
two "scattering angles", see section 4. In section 5 we derive a formula, suitable for numerical 
integration of the full matrix element, which we employ in the last section to demonstrate the 
convergence of the series forwards the exact result for a simple example. We conclude in section 7. 
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2 Reduction of the Collision Integral 



Omitting the superscripts denoting the particle species 1 in eqn. (4) we can write the collision inte- 
gral as: 

C[f](k) = 1 ^- [(2^5(E k +E p -E q -E r )5 3 (k+p-c l -r)\M\ 2 F[f] ]J * v , (5) 

k J v=p,q,r ^ > v 



where we introduced 



F[f] = (1 " r/fc)(l " efp)f q fr - fkfp(l - - Cfr) ■ (6) 



E v denotes the relativistic energy of the particles "v" on the mass shell, i.e. E v = \/v 2 + m 2 , 
with three-momentum v and mass m„. We write the 3 -dimensional <5-f unction as the Fourier 
transform of unity and switch to spherical coordinates: 



5 3 (k + p - q - r) = J eAC+P-i-r) ^ , (7) 



d 3 X 

( 

The collision term (5) then becomes 



(8) 



Here we have defined D as 



D(k,p,q,r) = ^ J dn p J dQ q J dQ r 5 3 (k + p-c l -r)\M\ 2 

= |^ J r X 2 d\ J e iXk dn x J e iXp dn p J e~ iAq dtt q j e~ iXr dtt r \M \ 2 . 

(9) 

Note, that this definition renders D(k,p,q,r) a dimensionless quantity. Due to the presence of 
the (5-function we expect that the result is non-zero only if q + r > \k — p\ and k + p > \q — r\, 
because the equation k + p = q + r does not have a solution otherwise, for whatever combination 
of the solid angles fl p , Q q , Cl r . Therefore, the result will be proportional to 2 

@(k,p,q,r) = @(q + r— \k — p\)@(k+p— \q — r\) 

= 0(min (k + p, q + r) — max (\k — p\ , \q — r|)) . (10) 



After computing D(k,p,q,r) we can proceed with the integration of the remaining energy 5- 
f unction in eqn. (8): 

c '[ / l (t) = 6sk/ /e^-^WW*,™,^^, <ii) 

'From here on we will always use the momenta k, p, q and r in connection with only one particle species, such 
that it serves as a label for the species at the same time. We also use the convention v = |v| if the distinction from the 
four-momentum is clear from the context. 

throughout we use the Heaviside-step-function 9 in the half-maximum convention (which will become relevant 
later). 
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where p = E£ — m| and E, p = E q + E r — E^. The Heaviside-f unction prevents us from inte- 
grating over combinations of q and r which are kinematically forbidden. Thus we have reduced 
the collision integral to a two-dimensional one, suitable for numerical integration. However, all the 
work is now hidden in the definition of D = D(k, p, q, r) which is characteristic for the scattering 
model, i.e. for the matrix element of the underlying theory for the scattering process under consid- 
eration. For completeness we present the analogous calculation for C 1 *^ 2 like collision integrals 
in appendix A. 

The computation of D is easily carried out for matrix elements squared with simple angular de- 
pendence such as the constant (|.M| 2 = const), for matrix elements in the Fermi approximation 
(\M.\ 2 oc (k • p)(q • r), (k • p), including re-namings of the momenta therein) and for resonant 
processes in the narrow width approximation (\M\ 2 oc S(s — m 2 x ), where mx is the mass of the 
particle in the intermediate state). However, in particle physics, one encounters matrix elements 
squared with a more complicated structure, such as products of tree-level s-, t- and u-channel 
contributions, for which the integrals D are in general unknown. 



3 A Simple Numerical Model 

In this section we establish a simple numerical model to benefit from the reduced form of the 
collision integral. For simplicity we assume a single particle species undergoing 2 — 2 scattering 
processes only. The system (1) then acquires the form 

L[f] = C[f] , (12) 

with C[f] from eqn. (11) and L[f] either from eqn. (3) or (2). 

In case the Liouville operator of the system is of the first kind with the extra term 

as compared to eqn. (2), which accounts for the expansion of the universe 3 , we first introduce the 
transformed variables 

x = Ma(t) , k = ka(t) , (14) 

with some suitable mass scale M and cosmic scale factor a(t). The Liouville operator in this new 
coordinates has the simpler form 

LI/]-*,M. (15) 

In the subsequent we omit the tilde over the transformed momenta and time, but it is important to 
remember that, in this case, the collision integrals need to be expressed in terms of the transformed 
variables as well. 



3 It is this term which prevents the Maxwell-Jiittner distribution function in general, from being an exact equilibrium 
solution for the Boltzmann equation in Robertson- Walker space-time. Only if the interaction rate is high enough its 
effect on the shape of the distribution function can be neglected. 
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Now we divide the physical relevant part of the positive real axis of momenta, i.e. we consider 
only momenta up to a maximum of k max , into a set of M disjoint (not necessarily equidistant) 
intervals Aki and choose a ki for each interval with ki G Aki (i = 1 . . . M). 
Then we make the approximation 

/ f(t, k) dk ~ f(t, hi) Ak t = UAh . (16) 

JAk t 

By integrating the left-hand side of (12) over the interval Aki we obtain 

Li= f ^-dk- y±Ah , or U = f Hx^ dk ~ Hx^Ah , (17) 
y Afc; eft dt J Akl dx dx 

For the right-hand side, we find 
1 M 



xDfap^kj) ^^ , (18) 



where p = yj (E ki + £; fe . - E h ) 2 - mj. 

This way, we turned the reduced Boltzmann equation into a coupled set of M ordinary differential 
equations for the discrete variables //: 

L l = C h (l = l...M). (19) 

Due to the finite momentum cutoff k max , this method is not conservative, i.e. energy and total 
particle number are not conserved inherently. In order to make the method energy and particle 
number conserving, the cutoff has to be chosen high enough. 

The number of equations, or grid points M, depends on the specific problem and the required 
accuracy. In any case eqn. (19) will represent a large system of differential equations. The amount 
of numerical work for the evaluation of the collision integral is of order 0(M 3 ). 
D can in principle be computed numerically, see section 5, and tabulated once and for all on 
the grid. For cosmological problems, however, the momenta remain to be scaled according to 
the time dependent connection (14), so that D needs to be re-evaluated in every time step. (A 
possible exception is the one of ultra-relativistic particles for which the scale invariance of D can 
be exploited.) This shows that the entire method relies on analytic expressions for D. In the 
following section a method is presented for the computation of D for arbitrary matrix elements. 



4 Expansion of the Scattering Kernel 

The spin averaged matrix element \Ai\ 2 will in general depend on Lorentz-invariant combinations 
of the four-momenta of the in- and outgoing particles, usually the Mandelstam variables s, t, u, 

s = (k+p) 2 , 

t = (k-q) 2 = m 2 k +m 2 q -2E k E q + 2\k\\q\cos(6 kq ) , 

u = (k - r) 2 = ra\ + m 2 r - 2E k E r + 2 |k| |r| cos (6> fcr ) . (20) 
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Figure 1: D 2 '°(2.0,p = qi + r — 2.0, qi, r) for qi < k. The flattening for r > k is common to all D 
All momenta are in relative units. 



ii.O 



In the following we take t and u as the independent variables and s is expressed by 4 



s = ^2 m ? — t — u . 

i=k,p,q,r 



(21) 



In order to evaluate eqn. (9) for arbitrary matrix elements we expand \Ai\ 2 in terms of cos (0kg) 
and cos (#fc r ) : 



fer ) 



(22) 



n=0 m=0 



Note that the coefficients A nm can depend on the magnitudes of the momenta. Upon integration 
of eqn. (9) we can then write 



D(k,p,q,r) = J]2 A nm {k,p,q,r)D nm {k,p,q,; 



(23) 



n=0 m=0 



assuming that the series converges for all relevant k, p, q and r (the momenta are still restricted by 
energy conservation). 



It can be advantageous to use different variables, in which case the results of this section can be adapted easily. 
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Figure 2: This plot continues fig. 1, D 2 '°(2.0,p = qi + r — 2.0, qi, r) for q l > k. 



In order to give a meaning to eqn. (23) we need to compute the integral 

D nm (k,p,q,r) = ^ J dn p J dQ q J dtl r 5 3 (k + p - q - r)(cos 9 kq ) n (cos 6 kr ) T 

X 2 d\ j e iXk dQ x J e iXp dn p 

J e - iAq (cos e kq ) n dn q J e - lXr (cos e kr ) m dn r . 



pqr 
64vr 5 



(24) 

Due to this definition the D nm, s are dimensionless, scale invariant functions, i.e. 
D nm (ak, ap, aq, ar) = D nm (k,p, q, r) for any a / 0. From the first line of eqn. (24) we can 
infer that, for given k, p, q, and r, the lowest order function D 0,0 (k,p,q,r) represents an upper 
bound for all higher order functions D nm {k, p,q,r). 

Before investigating further the general solution, we compute D nm for this simplest case, corre- 
sponding to \M.\ 2 = 1. 

We can evaluate all solid angle integrals, which \M. | 2 does not depend on, in eqn. (24), using 

47T 



e ±iXp dn p = ^ sm(Xp) 
Ap 

Thus D 0,0 simplifies to 5 

4 f°° dX 
D 0,0 (k,p,q,r) = - — / sin(AA;) sin(Ap) sinfAg) sin(Ar)— T . 

'■- 1 X z 



kir 



(25) 



(26) 



Without loss of generality we assume k, q, r > in the following. The cases k = 0, q = 0, r = can be 
understood in the limiting sense. 



7 



0.12 



0.08 



0.06 



0.04 



0.02 




qi = 0.2 



qio = 2.0 



20 



Figure 3: D 4 ' 3 (2.0,p = qi + r — 2.0, qi, r) for < k. The kink at r = k is common to all D nm (k, q 
r — k,q, r). For r < k — qwe have D nm = 0. All momenta are in relative units. 



Using the addition theorems for sin and cos, the result for this integral is found to be 

D°>°(k,p,q,r) = 

= — (\k — p + q + r\ — \k + p — q — r \ + \k + p + q — r \ + \k — p — q — r\ 

+ \k + p — q + r\ — \k — p — q + r\ — \k — p + q — r\ — \k + p + q + r\ \ , (27) 



or equivalently 



D '°(k,p,q,r) = — ( R{k - p + q + r) - R{k + p 



+ R{k+p + q — r) + R(k — p — q — r) 
+ R{k + p — q + r) — R{k — p — q + r) 
— R (k — p + q — r) — R (k + p + q + r)j , 



where we introduced the ramp function: 

, s , , fx, for x > 



(28) 



(29) 



Because of the proportionality of D nm to (k, p, q, r) from eqn. (10), we can infer that 

(k — p + q + r) > 0, (k + p — q + r) > 0, (k + p + q — r) >0 and (k — p — q — r) < for all 

values of k, p, q and r for which the prefactor is non-zero (none of the momenta can be greater than 
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Figure 4: This plot continues fig. 3, D 4 ' i (2.0,p = qt + r — 2.0, qi,r) for > k. 



the sum of the other three). Obviously, the term (k + p + q + r) is always positive. Introducing 
the abbreviations 



c\ = k + p — q — r, c 2 = k— p + q — r, c% = k — p 
and R 1 = R(c 1 ), R 2 = R{c 2 ), R 3 = R(c 3 ) , 



+ r 



for the remaining combinations with indefinite sign, we find for eqn. (28) the compact form: 

D°>°(k,p, q, r) = ^9 (k,p, q, r )(-R 1 -R 2 -R 3 + 2k). 



(30) 



(31) 



From this it is obvious that D°'°(k,p,q,r) < 1 everywhere. Since the smallest value, with all 
the Ri's being positive, is (k,p,q,r) (—k + p + q + r)/{2k) > it is also obvious that 
D 0,0 (k,p,q,r) > 0. Remembering the note from above, we conclude that \D nm (k,p, q, r)\ < 
D 0,0 (k,p,q,r) < 1 for all k, p, q and r. This important property guarantees pointwise abso- 
lute convergence of the series eqn. (23) if the series of the coefficients in the expansion (22), 
YlT=o Em=o A nm (k,p, q, r) , is pointwise absolute convergent. 

In case of massless particles eqn. (31) can be simplified further to (using energy conservation) 

D°'°{k,p,q,r) = -@(k,p,q,r)(k- R(q-k)- R(r-k)) 



k 



2k 



(k,p,q,r) (q + r — \q — k\ — \r — k\) . 



(32) 



We now turn to the computation of D nm for general n and m. In eqn. (24) the Fourier integrals, 



e ±iAq (cos 9 kq ) n dn q = I e ±iAq (kq) n dVL q , 



(33) 
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Figure 5: D 4 ' 3 (2.0,p = q + r — 2.0, q, r). fig. 3 and fig. 4 correspond to cuts with q = const (the thick 
line corresponds to qg = 1.8). 



on the unit sphere, can be expressed as a finite series of spherical Bessel functions of the first kind 
(see appendix B for the derivation): 

L-J 

J er^fa)" A*, = 4vr(±z)« E «n, a ^^(kA)- 2 ° , (34) 

with numeric coefficients 

(-l) a n! 

0-n,a — 71 77. TT • (35) 

a\[n — lay. 
Inserting eqn. (34) into eqn. (24) we find: 



- Lf J Lf J 

D nm = PV_j dXX 2 Y^(-i) n + m a n , a a mi(3 (2qX)~ a (2r\y 



i I »n.n"iii. H-'l") \ — ' ^ 1 / 

a=0 /3=0 



(36) 
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where the inner integral is again of type (34), such that we arrive at: 

[«J [SiJ ["+ m _( Q ,+ / 3)j 

D nm = ^fi2i2 2 £ (- 1 ) a+/3 «™,^ m ,/3«n +m -2( Q+ ^),K2g)- a (2r)-' 3 (2fc)- i 

a=0 0=0 Z=0 

x I(a + (3 + l;n + m — 2(a + /?) — /, 0, n — a, m — (3; k,p, q, r) , 

(37) 

with 

I(n;h,l 2 ,h,k;k,p,q,r) = / X 2 ~ n j h (kX)j l2 (pX)ji 3 (qX)ji 4 (rX)dX. (38) 

■/ o 

Unfortunately, the remaining integral over four spherical Bessel functions is known to represent 
a mathematical problem itself. Because of the rapidly oscillating integrand it is also diffcult to 
access by numerical methods. 

Due to Rayleigh's formula (A. 8), the integrand can always be decomposed into products of four 
sine and cosine functions and an inverse power of A, 



f 

Jo 



trig! (fcA) trig 2 (pA) trig 3 (qX) trig 4 (rA) ^ 
A™ 



where trigj (xX) is either sin(xA) or cos(xA). However the integrand in eqn. (39) has a non- 
integrable singularity at A = if the number of sines exceeds n. In principle, the problem can 
be circumvented by performing a Laurent series expansion of the integrand in eqn. (39) and sub- 
tracting the divergent part. The finite part can then be computed for all possible combinations of 
the indices. Since we expect a finite overall result for eqn. (37), the different divergent parts in the 
sum need to cancel. 

Here we make a different approach which is based on an explicit expression for the integral 

poo 

l(0;h,l 2 ,h,k;k,p,q,r) = / X 2 j h (kX)ji 2 (pX)ji 3 (qX)j k (rX)dX, (40) 

Jo 

valid, provided that there exists an integer number L which satisfies the conditions: 

\h ~ h\ < L < h + l 2 A |/ 3 - h\ < L < l 3 + h , and 

h + l 2 + L and Z 3 + Z 4 + L even . (41) 



In order to bring the integrals I{a + (3 + l;n + m — 2(a + (5) — 1, 0, n — a, m — (3; k, p, q, r) into 
the form of eqn. (40) we apply the recurrence relation for spherical Bessel functions [19]: 

3n(z) = 2n+1 tin-l(z) + j n +l(z)) ■ (42) 

Applying this relation r times with respect to j n (xX) yields a sequence of spherical Bessel func- 
tions of order n — r, n — r + 2 . . .n + r — 2, n + r and an overall prefactor (xX) r '. There- 
fore, we apply eqn. (42) I times with respect to j n+m _ 2 (a+(3)-l{kX), a times with respect to 
jn-a(qX) and (3 times with respect to j m _p(rX) in eqn. (36). This leads to a set of integrals 
of type 1(0; h,l 2 , 13, U; k,p, q, r), where < n + m — 2(a + (3 + I) < l\ < n + m — 2(a + (3), 
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h = 0, < n — 2a < I3 < n and < m — 2(3 < k < m. These integrals can be evaluated 
according to [20] if the conditions (41) on the indices are met. 

Since different authors found different, relevant expressions for integrals involving three Bessel 
functions [21-23], we repeat here the derivation for integrals involving four spherical Bessel func- 
tions from the former ones. 

With help of the closure relation for spherical Bessel functions (A.9), integrals of the form (40) 
can be reduced to integrals of three spherical Bessel functions: 

2 f'OO /'OO 

l(0;h,l 2 ,h,k;k,p,q,r) = - dXX 2 z 2 j h (kz)j L (Xz)ji 2 (pz)dz 

n Jo Jo 

poo 

x / z' 2 j h {qz')j L {Xz')j h {rz')dz' 
Jo 

= -/ dXX 2 I(h,L,l 2 ;k,X,p)I(h,L,k;q,X,r), (43) 

defining 

POO 

I(h,h,k;k,p,q) = / X 2 j h (kX)j h (pX)j h (qX)dX. (44) 
J 

Inserting the expression (A. 11) for l 2 ,k; k,p, q) found by Mehrem et al. [20] and performing 
the integration (after inserting the explicit representation (A. 12) for the Legendre polynomials) 
yields 

l(0;h,l 2 ,h,k;k,p,q,r) = 

(-i)-_ g ^V^TiP^e( fc> p >fl> r)( F ) (?) 

h h a -1 (h k A _1 v A ((2i 2 \ (2k\y 

0/ \ 0/ ^ 4^ \\ 2n J \ 2n 'JJ 

h l 2 — n l\ fh k — n ' I 
j l 



n=0n'=0 

h+h—n h+h—n' /, , /, , '7/ 

x E E i2/ ■ n;2/ ' • 11 d 



/=|; 1 -/ 2 +ri| l'=\l 3 -l 4 +ri 

L n l\(L n' l'\ (h l 2 L \ (h k L 
Oj \0 0j\n I l 2 -n)\n' I' k - n' 

J(k,p,q,r;n,n',l,l') 



k n q n ' 

where 



Jl 32 J3 

mi m 2 777,3 



(45) 



)-{£ 


J2 


J3l 




J5 





(46) 

denote the Wigner 3j- and 6j -symbols, respectively (these are purely numeric factors related to the 
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Clebsch-Gordan coefficients and Racah's W-coefficients, respectively [24]) and with 

J(k,p,q,r;n,n',l,l') = 

LfJ L^-J \-isi'-2t 

= EE a^^-w^' e e o - M c - 

s=0 t=o n=0 u=0 

x (k 2 - p 2 ) l - 2s -^(q 2 - r 2 ) l '- 2t ^U n+n , +2ll+ 2v + 2s+2 M '+i (k,p, q, r) , 

(47) 

where we have defined 

U a (k, p, q, r) = minCfc+p.g + rr-maxdp-fcMr-gir (4g) 

a 

The expected prefactor 6 (k,p, q, r) in eqn. (45) stems from the integration of the Heaviside step 
function in eqn. (A. 11). 

From eqn. (A.l 1) the result (45) also inherits the restrictions on the indices of the Bessel functions 
\h — fa\ < L < l\ + 12 A I/3 — Z4I < L < I3 + 14. Note, that, in general, eqn. (45) can be evaluated 
for different values of L and with different mapping of the indices l\, I2, h an d h, leading to 
different equivalent results. 

All sums in eqn. (37), eqn. (45) and eqn. (47) are finite, so they can be used to determine the 
functions D nm . In order to demonstrate the usefulness of this result we apply it explicitly to the 
cases of D°'° and D 2 ' . Evaluating (37) for n = m = we find 

D°'° = ^/(0;0,0,0,0;fc,p, g ,r). (49) 

TT 

Applying eqn. (45) leads immediately to 

6 (k,p,q,r) U x (k,q,p,r) 

2k 

Q(k,p,q,r) (min (p + r, k + q) - max (\-q + k\ ,\—r + p\)) 



2k 



(50) 



Differentiating the eight cases with sgn(q) = ±1, this can be shown to be equivalent to eqn. (31). 
A more sophisticated example is D 2 ' . Again, from eqn. (37), we find 



D 2 >° = ^(±WWAk,p,q,T)- 

TT \2 

/(-l;l,0,2,0;fc,p, g ,r) /(-!; 0, 0, 1, 0; k,p, q, r) \ 

2k + 2q J' K } 

We apply the recurrence relation (42) To the second and third term with respect to ji(k\) and 
ji(q\) respectively. This leads to 

D 2,o = 8p^r / 1 /(o _ ^ ^ ^ Q _ k ^ 9> r) + 1 J(0; 2) 0) 2> . fcjp> 9> _ (52) 
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The terms involving 1(0; 0, 0, 2, 0; k, p, q, r) did cancel exactly and both of the remaining integrals 
can be evaluated according to eqn. (45) (with the unique choice of L = and L = 2 for the first 
and the second integral, respectively), giving after some algebra: 

D 2,o Q(fc,P,g,r) 
8q 2 k 3 

x ({k 2 + q 2 ) 2 m (k, q, p, r) - 2 (k 2 + q 2 ) U 3 (k, q,p, r) + U 5 (k, q, p, r)) . 

(53) 



Eqn. (53) as well as all other functions D nm (k, p, q, r) can be brought into the compact form 

D^(k, p, q, r) = A^0- t (B lRl + B 2 R 2 + B 3 R 3 + C) . (54) 

We computed the numeric prefactor A and the momentum dependent coefficients, B\(k,p, q, r), 
B 2 {k,p, q, r), B 3 (k,p, q, r) and C(k,p, q, r) for n + m < 16. They are listed in appendix C for 
jjnm w j t j 1 n _)_ m < 5 6 coefficients B, L and C are homogeneous multivariate polynomials of 
degree 2 [n + m) and 2 (n + m) + 1 in fc, p, q, r. The number of elementary operations, necessary 
to evaluate D nm , increases with increasing n and m, however their shape permits considerable 
optimisation, especially when many D nm, s are to be computed. In addition, when dealing with 
networks of Boltzmann equations, they can be used for all matrix elements in the system. 

Fig. 1 and fig. 2 show D 2 '°(2.0,p = qi + r — 2.0, qi,r) plotted against r (all momenta are in 
relative units) for various values of qi. For r > k the graph of D 2 ' becomes constant. This can 
be understood by the observation that U a (k,q + r — k, q, r) is independent of r for r > k. The 
same relation holds for all other Z? n,0 's (and a corresponding one for D°' m ). Fig. 3 and fig. 4 show 
similar plots for D 4 ' 3 which do not have this property. Fig. 5 shows an surface plot of D 4,3 for 
fixed k, as a function of q and r, D 4,3 (2.0,p = q + r — 2.0, q, r). 

In general the shape of the graphs varies strongly with varying indices n and m. All functions 
possess a kink at r = k, because 

D^(k, P = q + r-k,q,r)= A *^±^ m (2 B 2 R(k - r) + 2 B 3 R(k -q) + C), (55) 

and the properties D nm (k,p, q, r) = for r < k — q, (0 (k,p, q, r) = @(q + r — k) = 
Q(p)) and lim D nm (k,p, q,r) = for k, p, and q held constant. The later one can be inferred 

from \D nm (k,p, q, r)\ < D°'°(k,p,q,r) and lim D 0,0 (k,p,q,r) = 0, which is obvious from 

r— »0+ 

eqn. (50). Similar relations hold for the q dependence (with k and r fixed) and in the case of 
massive particles. 



5 Numerical Integration of D 

In this section we derive a formula for D(k,p, q, r), suitable for numerical integration of arbitrary 
matrix elements. Again, we assume that the angular dependence of \M.\ 2 is given in terms of 

6 They become lengthy for greater indices. 
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COS 

s 



ds (0k q ) and cos (9k r ), i.e. in terms of the momentum transfer t and u. A possible dependence on 
s can be expressed in terms of t and u exploiting energy and momentum conservation, eqn. (21). 

We orientate the coordinate system such that the z-axis points in the direction of k. We can then 
write 

cos (#Av) = ^ " v = cos6*acos6* v + smQ\sva.Q v cos{(j)\ — <p v ) , 

cos (6 kq ) = COS 6> q , 

cos (9kr) = cos9 r . (56) 
Using eqn. (34) for the $7 P integration we can write eqn. (9) as 

D(k,p,q,r) = JgLj d*\ e ^^_M J e -iM dnq J e -^ d n r \M\ 2 . (57) 

We can then perform the integration over <fi g and 4> r since \M | 2 does not depend on these angles 

J e-' lXr dn r \M\ 2 = 

e -i\r cos 6 X cos 6 r dcos Q r J e -iArsine A sine r cos(0 A -<M j^) 2 ^ 

= [ e~ iXrcose * cosd r\M\ 2 dcos9 r f 2cos (Ar sin A sin r cos (cp r )) dcf) r 
J-i Jo 

= 2vr ^ e- iXrcosd * cosd rJ (\rsm9xsm9 r ) \M\ 2 dcos9 r , (58) 

where we have taken into account the fact that the inner integral does not depend on 4>\, since the 
integration is over 2ir and the odd symmetry of the imaginary part of this integral with respect to 
cf> r . The remaining integral was recognized as the integral definition of Jo, the Bessel function of 
the first kind of order zero, see e.g. [19]. 

Inserting eqn. (58) into eqn. (57) gives 

D»,M, T )-%£°\'*iMldX, (59, 

with 

I = J e iXk cos dx dn x J 1 e- iXq cos e * cos e « d cos # q J 1 e- iXr cos e * cos dr d cos 9 r 
x J (\qsm6\siiii9q i ) Jo(\r sin 6>a sin 6* r ) \ M\ 2 . 

(60) 

For the product of the two Bessel functions we may use the relation (A.7). Inserting it into eqn. (60) 
and interchanging the order of integration, we find 



= 2^y J d cos 9 n d cos 9 r \M\ 2 I' , (61) 
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with 

jt _ I I ^i\(k— gcos # q — r cos 6 r ) cos 0\ 

n Jo Jo 



x Jo ^A sin 9 x y (q sin # q ) 2 + (r sin 9 r ) 2 — 2qr sin 6> q sin 9 r cos (x) J sin 9\ dx d9\ . 

(62) 

Again interchanging the order of the integrals and exploiting the odd symmetry of the imaginary 
part of the exponential, we get 

1 f 7 ^ f 7 ^ 

I' = — / cos (A (k — qcosOq — r cos6> r ) cos6>a) 
* Jo Jo 

x Jo ^A sin 9\\J (q sin # q ) 2 + (r sin 9 r ) 2 — 2qr sin # q sin 9 r cos (x)^j sin 6*^ d9\ dx , 

(63) 

Now we apply the integral (A.6) to arrive at 



2 r sin ( A v7M) 

I' = - v - ' dx , (64) 

vr 7o AyTR 

with 

/(cos x) = (k — q cos # q — r cos 6* r ) 2 + (g sin 6* q ) 2 + (r sin # r ) 2 — 2 qr sin 6* q sin # r cos x , (65) 

considering < 9 q , 9 r < ir as parameters. If we interpret x as the polar angle enclosed by q and 
r we can write, 

/(cos a;) = (k - q) 2 + (k - r) 2 - (q - r) 2 - k 2 + q 2 + r 2 

= (k - q) 2 + (k - r) 2 - (q - r) 2 - k 2 + q 2 + r 2 + (£ fc + E p -E q - E T f 
= s + t + u-^m 2 +p 2 , (66) 

under the assumption of energy and momentum conservation. 

Since qr sin 6> q sin 9 r > for all q , 9 r , the function f(x) takes its minimum for cos x = 1: 

/(l) = (A; - q cos 6> q - r cos 9 r ) 2 + (q sin (9 q - r sin (9 r ) 2 > . (67) 
I' is therefore well defined. 

Inserting both, (61) and (64), into eqn. (59) and again exchanging the order of integration, we find: 

•1 r l 



pqr f f 

D(k,p,q,r) = — ^ / / dcos9 cl dcos9 r 
71 J-iJ-i 



roo gm , 



/ \„\ sin f A \/ f (cos x) 
x \ M f I dx I X 2 ^tM \ V_ ^d\. (68) 
lo Jo Ap X^ficosx) 
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In the rightmost integral we recognize the closure relation for spherical Bessel functions (A.9) for 
n = 0. 
This leads to 

D(k,p,q,r) = ^- [ [ dcosO^dcos 9 r \M\ 2 1" , (69) 



where we defined 



(70) 



(71) 



with 

F(e q ,e r ) = (p?-/(i))(/(-i)_p=«) 

= (2qr sin # q sin 9 r ) 2 — {k — q cos # q — r cos 6> r ) 2 + ((? sin # q ) 2 + (r sin 6* r ) 2 — 

Because of eqn. (66) the 5-f unction in eqn. (70) ensures energy and momentum conservation. 
The final expression for D reads: 

D(k, p,q,r) = ^ [ j^ 1 d cos fl q d cos 9 r , (72) 
^ JA \J r [Uq, U r ) 

where the domain of integration A is given by — 1 < cos 9 q , cos 9 r < 1 and F(6 q , 9 r ) > 0. Note, 
that the term 1 / y / F(9 q , 9 r ) seems to be absent in an analogous expression in [1]. We are confident 
that it must be there, because in the simple case of = 1 the result would be iqr/n, without 
F, which is different from eqn. (31). 

The expression (72) for D has only a two-dimensional integral and is by far superior to eqn. (57) 
with respect to numerical integration. The integrand may have singular points at the boundary of 
A. Hence the routines for numerical integration must be chosen adequately. Usually, for a numer- 



ical method, it is sufficient to know D(k,p = y(E q + E r — E^) 2 — m%, q, r) for a finite set of 

momenta {ki, qj, r{\ on a grid. Therefore it is possible, in principle, to tabulate D through numer- 
ical integration of eqn. (72). As pointed out above, for applications in cosmology, this relation is 
only of restricted use, since the momenta or, equivalently, the particle masses are scaled in each 
step of the time evolution, so that the values of D(ki,p, qj,ri) need to be recomputed permanently. 



6 Convergence of the Method 

To demonstrate the application and convergence of the method described above, we apply it to a 
(hypothetical) matrix element involving tree-level t- and u-channel contributions. We compare the 
exact numerical result (by the formula derived in section 5) with the approximate analytical result 
according to the truncated series expansion of section 4. 

Without specifying a theory we take the matrix element to be (for simplicity we assume that all 
in- and outgoing particles are massless): 

\M\ 2 =( 9 +^ V=f- " + ~ " V, (73) 

\mx 2 — t nix 2 — u J \ 2 kq (a — cos q ) 2 kr (b — cos 9 r ) J 
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Figure 6: D(k, q + r — k,q,r) with mx = 20.0, k = 15.0, q = 10.0. Exact numerical result D num , result 
in large mx limit -Dfermi, successive analytic approximations Dq ■ ■ ■ Dg (D n+m corresponds to 



an expansion of | M. | up to order n + to in the cosines). 

= 1 + m x 2 



with a = 1 + mx 2 / {2 kq) > 1 and 6 = 1 + mx 2 / (2 kr) > 1 and some mass mx 
intermediate state and a coupling g with mass dimension 2. 

A Taylor series expansion in cos 6> q and cos r leads to: 

W 2cos^ V_£_I( 1 + 2c««, >| 

1 / cos 8„ cos#- (cos 9^) 



of the 



Q 2 1 



2 /c 2 qr afo 



9 ■ + _ 

a 



COS# r (cOS0 q ) 2 COS 6> r COS 6> q (COS 

6 a 2 a& 6 2 



. /•) is found by 



The angle-integrated matrix element D(k,p,q, 
(cos# q ) n (cos# r ) m by D nm : 

[k 2 r 2 b 2 \ b 

J 2 ' z> 



(cosfl r ) 2 \ 

(74) 
of 



+ 



9 



substituting every appearance 
1 / 2D 1 ' \ g 2 1 / 2D ' 1 \ 

^{ 1 + — + --j + u^v 2 { 1 + — + --j + 

1 / D 1 ' D ' 1 D 2 ' D 1 ' 1 D°< 2 \ 



\M\ 2 = 

Ak 2 q 2 a z \ a 

,2 i / D l,0 D 0,l ^2,0 D l,l jy 

b a 2 ab b 2 



2 k 2 qr ab 



Setting cos # q = 1 and cos 



r = 1 yields the series of coefficients which convergt 

cose q =cos6» r =i \2 kq (a — 1 



;es to 



g 1 
2kr (6-1), 



(75) 



(76) 
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Figure 7: D(k, q + r — k,q,r) with mx = 10.0, k = 10.0, q = 15.0. Exact numerical result D nam , result 
in large mx limit Df erm \, successive analytic approximations Dq, . . . D$. 



According to what has been said above \Ai\ D°'°(k,p,q,r) represents an upper 

COS 0q=COS 9 r = l 

bound to D(k,p, q, r). On the other hand, in the low energy limit, eqn. (73) can be expanded in 
terms of inverse powers of mx 2 '- 



\MY 



9 



m 



x 



4 + 4- 



ii 



+ 5 



m x 



+ 2 



tu 

m x L 



+ 5- 



(77) 



which corresponds to the Fermi-approximation and includes only (cos 6> q ) n (cos 9 r ) m terms with 
n + m < 2 at this order (i.e. these terms can be integrated as in the previous literature). 

fig. 6 and fig. 7 show the graphs of the exact numerical result (eqn. (72)), the graphs for the "Fermi- 
approximation" (eqn. (77)), the graphs of the theoretical upper bound (eqn. (76)) and the graphs 
corresponding to successive approximations according to eqn. (75) for two sets of parameters. 
Fig. 6 and fig. 7 show that for parameters, for which the expansion in inverse powers of mx 2 
naturally fails, the approximation by the truncated D nm expansion gives quite good results. 

The rate of convergence of successive ZVs torwards the exact result depends on the momenta 
since the coefficients in the expansion (74) are momentum dependent. In this case it becomes 
worse for mx 2 / (2 kq) <C 1 and mx 2 / (2 kr) <C 1. 

7 Conclusion 



In state-of-the-art computations in astroparticle physics usually all species apart from the neu- 
trinos, which experience only weak interactions, are assumed to be in exact kinetic equilibrium 
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and heterogeneous (at best) networks of Boltzmann- and rate equations are solved instead of the 
full system of kinetic equations. The number of species involved in realistic systems is large and 
requires a unified treatment of the different particles and interactions. 

A method for the solution of the space homogeneous Boltzmann equation (with isotropic distribu- 
tion function), for general scattering laws was presented here. The method relies on the expansion 
of the matrix element in terms of cosines of two "scattering angles". For the separate terms in this 
expansion the full angular integration was carried out. The functions D°>°, D ' 1 and D ' 2 , corre- 
sponding to matrix elements in Fermi-Approximation, were used in previous literature to compute 
non-equilibrium corrections to the neutrino-distribution functions and have been obtained in the 
lowest order of the expansion. 

Though our starting point was the relativistic form of the Boltzmann equation, as encountered in 
astroparticle physics, the method can be used for the non-relativistic equation as well. In any case, 
it allows for the full angular integration of the scattering kernel, reducing the collision integral 
from effective dimension 5 to dimension 2. The only prerequisite is that the matrix element can be 
expanded into a series of the scattering angles and that this series converges rapidly enough. The 
quantum statistical terms for blocking and stimulated emission can be carried along. 

In the introduction we mentioned that, at high densities and temperatures, modifications of the 
Boltzmann equation might become necessary (if these are sufficient at all). One such modifica- 
tion, which has been suggested on various occasions, is the inclusion of higher order scattering 
processes. Since the representation of the angular integral in terms of spherical Bessel functions 
(36-37) and the integral (45), by means of eqn. (43), can be generalized for higher order processes 
(in which case more than four Bessel functions appear in the integrals) the method, described 
above, can in principle be used to reduce the corresponding collision integrals from dimension 
3 (n — 1) to n — 2. (n is the number of particles involved) 

The functions D nrn , though of very simple structure, can become lengthy for higher orders. More- 
over, due to the presence of very different relaxation time-scales the system of ODE's, correspond- 
ing to the numerical method presented in section 3, tends to behave stiff. Therefore, in possible 
implementations, careful optimisation for efficiency and stability is necessary. 

Let us add, that the expansion of the scattering kernel in terms of the cosines of the angles is not 
a new idea. For example in [25] an expansion of the scattering kernel has been combined with a 
moment method for the non-relativistic, inhomogeneous Boltzmann equation. The expansion of 
generic kernels with full integration of the angular part, in the space-homogeneous and isotropic 
case, seems to be new, however. 
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A Reduction of C 1 ^ 2 like Collision Integrals 

For collision integrals describing decays and inverse decays the angle integrated matrix element 
can be defined similarly to eqn. (24). In this case the matrix element is a constant and only the 
zeroth-order integral, corresponding to \M\ 2 = 1, has to be computed. The collision integral 
reads 

C^lJ] = ^ J(2*)W> (k - , - ,) \M\ *FV\ (2n % t (2n f 2Er . (A.) 

with F'{f] = (1 - Cfk)f q fr ~ /fc(l - - C/r)- Performing the same steps as in the main 

text, we derive 

where E q = Ej, — E r , q = \J E 2 — m 2 and we have defined the function D' as 

D'(k,q,r) = J ' X 2 dX J e iXk dQ x J e~ iXci dn q J e- iXr dQ r \M\ 2 . (A3) 



For \M\ 2 = 1 we find 



8 /*°° 

D'(k,q,r) = — I sin(AA;) sin(A(7) sin(Ar) 
™ J 



dX 
T 



= — — (sgn(r — q — k) — sgn(r — q + k) — sgn(r + q — k) + 1) . (A.4) 
k 

B Relations Involving Bessel Functions of Integer and Fractional Or- 
der 

For the reader's convenience we collect some facts about Bessel functions of the first kind J n and 
spherical Bessel functions of the first kind j n . They are related by 

3n{z) = y^4+i/2(2) • (A.5) 

In the main text we employ the following integral of Bessel functions from [26] : 

fr sin ( y/a 2 + (3 2 ) 

/ cob {P cob 6) J {a wad) an Ode = 2 V ' , (A.6) 

Jo v a 2 + (3 2 

and for the product of two Bessel functions: 

J (a) J (/?) = ^ J J (Va 2 + /3 2 -2a/3cos (x)) . (A.7) 

Using Rayleigh's formula, the spherical Bessel functions can be computed iteratively from the 
sine function: 

= *"(-;*)" ^ 
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They satisfy the closure relation [27] 

v 2 poo 
'0 



2z 2 f'°° 

— / A 2 j„ (Xz) j n (Xz') dX = 5(z- z') . (A.9) 
71 Jo 



Several authors have derived expressions or algorithms for the computation of integrals involving 
products of three spherical Bessel functions [21-23] 

poo 

I(h,l 2 ,l 3 ;k,p,q) = / X 2 j h (kX)j h {pX)j h (qX)dX. (A.10) 
J o 

Here we cite the explicit result, found by Mehrem et al. [20] by relating I(li,l 2 ,h',k,p,q) to 
known integrals of three spherical harmonics: 

I(h,l 2 ,l 3 ;k,p,q) = 

K®( q -\k-p\)Q{k+p-q)^- 1 * / —— f h h Z 3 \ 1 

x^(2/ 3 2n) 1/2 ( P Ar E ( 2/ + 1 ) 

n=0 l=\h-(.ls-n)\ 

h h-n l\ (h n l\ fh h h\ p ( k 2 + p 2 - q 2 



0J \0 0J [n l 3 -n I j 1 \ 2kp 

(A. 11) 



where Pi denotes the Legendre polynomials with explicit representation [19] 

2H\{l-i)\{l-2i)\ 



Pi(x)-T ( " 1)t(2j " 2i)! (A 12) 



i=0 

Wigner's 3j symbol is related to the Clebsch-Gordan coefficients by [24] 



mi m 2 m 3 J ^2j 3 + 1 



-(jimij 2 m 2 \j 3 -m 3 ) (A.13) 



Wigner's 3j symbols are equal to zero, unless mi + m 2 + m 3 = 0, \rrii\ < y\ and \ j\ — j 2 \ < j 3 < 

ji + 32 

Wigner's 6j symbols are related to Racah's W-coefficients by 

\ jl J - h ) = (-iy i+32+u+35 W(jd23di;hk). (A.14) 

[ 34 35 36 J 

In the main text, eqn. (34), we use an expansion into spherical bessel functions [28,29]. We start 
with the expression (given in [29]): 
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where x, fj G R d with (fj\fj) = 1, j u (z) = T(u + 1) (|)^ J v (z) and the integration is over the 
(d — l)-sphere S"* -1 . P n is a homogeneous polynomial of degree n on R d and (.|, ) denotes the 
inner product. For d = 3, we find 



/ 



I - 1 

/ j\ n 2 (-'\) a /9 \ 

e^W) dfy = J E (x J 3n-a(x)(A a P n )(x) . (A.16) 



o=0 

Choosing P ri (r ? ) = (kr)) n we get with (A Q P n )(x) = n!/(n - 2a)! • (kx) n " 2a 

LfJ 

a! (n - 2a)! (2x) 
Substituting x — > ±Ar/ reproduces eqn. (34). 



y e ix ^(kf/) r 



C The Integrals D nm 

The functions D nm (k, p, q, r) can all be written in the form 

D^(k, p, q, r) = A ^l^l (B lRl + B 2 R 2 + B 3 R 3 + C) . (B-l) 

In the following we list the coefficients A, B\ and C, which themselves depend on the momenta, 
for D nm with n + m < 5 and n < m (D n ' m with n > m can be derived from £) m > n by 
interchanging q and r). The expressions B 2 (-B3) are found by substituting in Pi the term c\ by 
c 2 (c 3 ) and fi by « (J;™) for alH. 

D°>°(k,p,q,r) : (B-2) 
4 = 1/2, C = 2fc, 
Si = -1, 

D°>\k,p,q,r) : (B-3) 
4= -1/12, C = -4fc 3 , 
Si = /2C1 - ci 2 + /1 , 

[/1 = 6kr,f 2 = 3fc-3r] 

D°> 2 (k,p,q,r) : (B-4) 

A = iio ' c *= 8fc3 ( 2fc2 + 5r2 )> 

-Bi = / 2 ci + /3C1 2 + / 4 ci a - 3 ci 4 + /1 , 

[ /1 = -60fcV,/ 2 = -60kr(k-r),f 3 = -20 r 2 - 20 fc 2 + 60 fcr, / 4 = -15r + 15fc] 
D°' 3 (k,p,q,r) : (B-5) 
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A = ~560' C=-16fe 5 (2fc 2 + 7r 2 ), 

Bi = feci + / 3 ci 2 + / 4 ci 3 + fed 4 + fed" - 5ci 6 + fe . 



[ fe = 280r 3 fc 3 ,/ 2 = 420fcV(fc -r),/ 3 = 140 kr(2k — r) (k - 2 r) , fe = 70 (fc - r) 
x(r 2 - 5fcr + fc 2 ),/ 5 = -42 (2fc - r) (k-2r),fe = -35 r + 35fc] 



D u >*(k,p,q,r) : (B-6) 
A = TMo' c = 32fc5 ( 36fcV + 63 r 4 + 8fc 4 ), 

Si = fed + fed 2 + / 4 ci 3 + fed 4 + feci 5 + fed 6 + feci 7 - 35 ci 8 + fe , 

[ fe = -5040fcV,/ 2 = - 10080 fcV (fc -r),fe = -3360fcV(3r 2 + 3fc 2 -7kr), 
fe = -2520 fcr (fc -r) (2r 2 - 7 fcr + 2 fc 2 ) , / 5 = 10080 fc 3 r - 1008 r 4 + 10080 fcr 3 - 1008 k 4 
-19656 fcV,/ 6 = 840 (k - r) (2r 2 - 7fcr + 2fc 2 ),/ 7 = 2520 fcr - 1080 fc 2 - 1080 r 2 , 
/s = -315 r + 315fc] 



D°> 5 (k,p,q,r) : (B-7) 
A= "44^52 ' C = ~ 64fe7 ( 44fe2r2 + 8fc4 +99r 4 ) , 

Bi = /2C1 + / 3 Ci 2 + feci 3 + feci 4 + feci 5 + feci 6 + feci 7 + feci 8 + feoci 9 - 63ci 10 + fe , 

[ fe = 22176 r 5 fc 5 ,/ 2 = 55440 fcV (fc - r) , fe = 18480 fc 3 r 3 (4r 2 - 9 kr + 4 fc 2 ) , fe = 55440 fcV 
x(fc - r) (r 2 - 3fcr + fc 2 ),/ 5 = 11088 fcr( - 14fc 3 r + 2r 4 - 14fcr 3 + 25fcV + 2fc 4 ), 
fe = 1848 (fc - r) (2 r 4 - 28 fcr 3 + 67 fcV - 28 fc 3 r + 2 fc 4 ) , / 7 = -99000 fc 2 r 2 - 7920 fc 4 
-7920 r 4 + 55440 fcr 3 + 55440 fc 3 r, fe = 6930 (fc - r) (r 2 - 3 fcr + fc 2 ) , fe = -3080 fc 2 - 3080 r 2 
+6930 fcr, /10 = -693 r + 693 fc] 



D l ^(k,p,q,r) : (B-8) 

^ = 120' c ' = 4fe3 ( 3fc2 + 5 P 2 - 5 9 2 - 5r2 )' 
Bi = /2C1 + / 3 ci 2 + feci 3 - a 4 + fe , 

[ fe = -60fcV,/ 2 = -30fc( - 2gr + fcg + fcr),/ 3 = 20 fcg + 20 fcr - 20 gr - 10 fc 2 , 
/ 4 = -5g-5r + 5fc] 



D^{k,p,q,r) : (B-9) 

A = > C = -16fc 5 ( - 14g 2 + 14p 2 + 7r 2 + 4fc 2 ) , 

1680 v ' 

Bi = feci + feci 2 + / 4 ci 3 + feci 4 + feci 5 - 3 ci 6 + fe , 



[ fe = 840gr 2 fc 3 ,/ 2 = 420fc 2 r(fcr - 3gr + 2fcg),/ 3 = 140 fc(2 rfc 2 + 6 r 2 q - 3 r 2 fc + 2 gfc 2 - 8rgfc), 
-56fc 2 + 126fcg,/ 6 = 21fc-21r-21g] 



/ 4 = -280 gfc 2 + 630 rgfc + 70 fc 3 - 280 rfc 2 - 210 r 2 g + 210 r 2 fc, / 5 = 126 fcr - 126 gr - 42 r 2 



D l > 3 (k,p,q,r) : (B-10) 
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A = , C = 16fc 5 (54fc 2 p 2 -63gV - 54g 2 fc 2 - 63r 4 + 27fcV + 63pV + 10fc 4 ) , 

Bi = ha + / 3Cl 2 + / 4Cl 3 + / 5Cl 4 + / 6Cl 5 + / 7 ci 6 + fsc-J - 5 Cl 8 + /i , 

[ /i = -5040r 3 fc 4 g,/ 2 = -2520 fcV (fcr - 4 gr + 3 fcg) , f 3 = -840 fc 2 r ( - 4 r 2 fc + 12 r 2 g + 3 rk 2 
-18rgfc + 6gfc 2 ),/ 4 = -1260 fc( - 4 r 3 g + 11 kqr 2 - 7fc 2 gr + gfc 3 - 3 fcV + fc 3 r + 2 kr 3 ) , 
/ B = 6048 fcgr 2 + 1764 fc 3 r + 1008 kr 3 + 1764 gfc 3 - 1008 r 3 g - 252 fc 4 - 6804fc 2 gr - 2772 fcV, 
/ 6 = 1008 r 2 k + 2520 rgfc - 1008r 2 g - 1134gfc 2 - 1134rfc 2 - 168r 3 + 294fc 3 ,/ 7 = 360 kr 
-360 gr- 144r 2 + 360 fcg- 162fc 2 ,/ 8 = 45 fc - 45g - 45r] 



D^(k,p,q,r) : (B-ll) 
A = — , C = -64fc 7 (99r 4 +88fc 2 r 2 + 396 p 2 r 2 - 396 g 2 r 2 + 176 k 2 p 2 - 176g 2 fc 2 + 24 fc 4 ) , 
Si = / 2 ci + / 3Cl 2 + / 4 ci 3 + / 5Cl 4 + / 6Cl 5 + / 7 ci 6 + / 8 ci 7 + / 9 ci 8 + /ioci 9 - 35ci 1() + /i , 

[ /i = 110880 r 4 gfc 5 ,/ 2 = 55440 fcV ( - 5 qr + kr + 4 kq) , f 3 = 18480 fcV (20 r 2 g - 5 r 2 k 
-32rgfc + 4rfc 2 + 12gfc 2 ),/ 4 = 18480 fc 2 r (5 kr 3 - 15 r 3 g - 8 fcV + 41 fcgr 2 -30fc 2 gr 
+3fc 3 r + 6gfc 3 ),/ 5 = 3696 k( - 15r 4 fc + 30r 4 g - 66fc 3 gr + 6gfc 4 - 141fcr 3 g + 6fc 4 r 
-30 fc 3 r 2 + 171 fc V 2 + 41 fcV) , f 6 = 18480 r 4 fc - 40656 fc 4 r + 240240 fc 3 gr + 105336 fc 3 r 2 
+ 184800 kr 3 q - 18480 r 4 g - 86856 fcV + 3696 fc 5 - 40656 gfc 4 - 382536 fc 2 gr 2 , f 7 = -2640 r 4 
-126720 fc 2 gr - 26400 r 3 g + 118800 fcgr 2 - 5808 fc 4 + 34320 k 3 r + 34320 gfc 3 + 26400 kr 3 
-54648 fcV,/ 8 = 14850 r 2 fc - 14850 r 2 g - 15840 rk 2 + 34650 rgfc - 15840 gfc 2 + 4290 fc 3 
-3300r 3 ,/ 9 = -1760 fc 2 +3850 fcg - 1650 r 2 - 3850 qr + 3850 kr, f w = -385 r - 385 g 
+385 fc] 



D 2 > 2 (k,p,q,r): (B-12) 

A = — — , C = 16fc 5 (28g 2 r 2 + 7r 4 + 22 fc 2 p 2 + 7p 4 + 2 g 2 fc 2 + 3 fc 4 - 14pV 
3360 

-14p 2 g 2 + 2fc 2 r 2 + 7g 4 ) , 
Bi = /2C1 + /3C1 2 + / 4 ci 3 + / 5 ci 4 + / 6 ci 5 + / 7 ci 6 + f 8 d 7 - ci 8 + /1 , 

[ /1 = -1680g 2 fc 4 r 2 ,/ 2 = -1680fc 3 gr( - 2 gr + fcg + kr), f 3 = -560fc 2 (fcr - 3gr + fcg) 
x ( - 2 qr + kq + kr) , / 4 = -140 fc(l2 gr - 5 kr - 5 fcg + 2 fc 2 ) ( - gr + kr + fcg) , 
/ B = -336 gV + 1008 fcgr 2 + 1008 fcg 2 r - 476 g 2 fc 2 + 336 fc 3 r - 1344 fc 2 gr + 336 gfc 3 
-476 fcV - 56fc 4 ,/ 6 = 56 ( - g - r + fc) (3gr - 3 fcg + fc 2 - 3 kr) , f 7 = -72 gr - 24 g 2 
-24 r 2 - 32 fc 2 + 72 fcg + 72 kr, fa = -9 g - 9 r + 9 fc] 



D 2 > 3 (k,p,q,r) : (B-13) 
A = - 1 , C = -32 fc 7 (297 p 4 + 198 pV - 22 g 2 fc 2 + 297 g 4 + 66 fc 2 r 2 + 462 fc 2 p 2 

-495 r 4 + 396gV + 41 fc 4 - 594 p 2 g 2 ) , 
Bi = faci + /3C1 2 + / 4Cl 3 + / 5Cl 4 + / 6Cl 5 + / 7 ci 6 + / 8 ci 7 + / 9 ci 8 + /ioci 9 - 15ci 1() + /1 , 
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[ /i = 110880 q 2 r 3 k 5 ,f 2 = 55440 k 4 qr 2 ( - 5 qr + 3 kq + 2 kr) , f 3 = 18480 k 3 r (20 q 2 r 2 + 6 k 2 qr 
+6 g 2 fc 2 - 21 kq 2 r - 12 fcgr 2 + 2 fcV) , / 4 = 27720 fc 2 (2 fc 3 gr + 17 kq 2 r 2 - 8 k 2 qr 2 - 10 g V 
+g 2 fc 3 + fc V + 9 fcr 3 g - 2 fcV - 8 k 2 q 2 r) , f 5 = 5544 k (9 fc V - 8 fc V + 20 g V + 45 k 2 qr 2 
+2k 4 r - 57 kq 2 r 2 - lSk 3 qr - 29kr 3 q + 2qk 4 - 8 g 2 fc 3 + 41 fc 2 g 2 r) , / 6 = 110880 fcgV 
-152460 k 2 qr 2 + 99792 fc 3 gr + 37884 g 2 fc 3 - 18480 g 2 r 3 - 130284 fc 2 g 2 r - 16632 k A r + 41580 fc 3 r 2 
-26796 fcV + 55440 fcr 3 g + 1848 k 5 - 16632 qk 4 , f 7 = -21780 fc 2 r 2 - 2376 k 4 - 7920 r 3 g 
+ 14256 gfc 3 - 53856 fc 2 gr + 39600 kq 2 r - 18612 g 2 fc 2 + 47520 fcgr 2 + 7920 kr 3 - 15840 gV 
+ 14256 k 3 r,f R = 198 ( - g - r + fe) ( - 25 kq + 25 qr + 9 k 2 - 25 kr + 5 r 2 ) , / 9 = -660 r 2 
-748 fc 2 - 550 g 2 + 1650 kr + 1650 kq - 1650 gr, /i = -165 g+ 165 fc - 165 r] 

The D nm, s do not possess any singularities inside the domain of integration, which is an almost 
essential feature for the numerical solution of the Boltzmann equation. Note, however, that the 
expressions given here are neither optimised for numerical efficiency nor for stability. Not all of 
the terms B{ need to be computed in every step since Ri = R(ci) = for q < 0, and quantities 
such as powers of the momenta, which appear several times, need to be computed only once for 
all D nm, s. In an implementation according to the model presented in section 3 the D nm, s need 
to be computed only once for all matrix elements in the system (with the possible exception of C 
and the q terms including p which is determined by energy conservation). 
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